%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 禁忌搜索算法
% 简介： 禁忌搜索算法寻找最大值
% 作者：Zhaojiang
% 日期：2023/10/10
% 企鹅：277746470
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% target fucntion
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
pause(2);
close all;clear;clc;
%% 设定算法参数
rng(1234);
var_num = 2;% 2个求解变量
lower_bound = -4;% 解的最小值（小数后4位）
upper_bound = 4;% 解的最大值（小数后4位）
bound_length = upper_bound-lower_bound;
omega = 1;%自适应权重
omega_decrease_rate = 0.9;% 权重衰减值(<1)
near_var_num = 5;% 邻域解生成个数
max_iteration = 200;% 最大迭代次数
eps = 1e-5;% 禁忌表精度
tabu_length = randi([near_var_num+1,2*near_var_num],1,1);% 禁忌表长度
%% 算法数据初始化
tubu = zeros(tabu_length,var_num);% 禁忌表初始化
now_tabu_length = 1;% 禁忌表当前禁忌表长度
var_value = lower_bound+bound_length*rand(1,var_num);% 当前变量值
% var_value = [0;2];
func_value = target(var_value(1),var_value(2));% 当前函数值
near_var_value = zeros(near_var_num,var_num);% 邻域变量初始化
near_func_value = zeros(1,near_var_num);% 邻域函数值初始化
best_near_var = var_value;% 邻域最佳变量
best_near_func_value = func_value;% 邻域最佳函数值
best_iteration = 0;% 最佳迭代次数
best_var= var_value;% 最佳变量值
best_func_value = func_value;% 最佳函数值
all_best_func_value =zeros(1,max_iteration);% 最佳适应值记录
in_tabu = zeros(1,var_num);% 禁忌真值表
%% 开始迭代求解
for iteration=1:max_iteration
    %自适应权重
    omega = omega * omega_decrease_rate;
    for i = 1:near_var_num 
        % 生成邻域解
        near_var_value(i,:) = var_value+(2*rand(size(var_value))-1)*omega*bound_length;
        % 限幅
        for j= 1:var_num
            if near_var_value(i,j) > upper_bound
                near_var_value(i,j) =upper_bound;
            end
            if near_var_value(i,j) < lower_bound
                near_var_value(i,j) =lower_bound;
            end
        end
        % 计算邻域函数值
        near_func_value(i) = target(near_var_value(i,1),near_var_value(i,2));
        % 是否为邻域最佳函数值
        if near_func_value(i) > best_near_func_value
           best_near_func_value = near_func_value(i);
           best_near_var = near_var_value(i,:);
        end  
    end
    % 历史最佳值、领域最佳值、当前值之差
    delta1 = best_near_func_value - func_value;
    delta2 = best_near_func_value - best_func_value;
    if delta1 <= 0 % 领域最佳值<当前值
        % 将领域最佳值写入禁忌表
        tubu(now_tabu_length,:) = best_near_func_value;
        now_tabu_length = now_tabu_length + 1;
    elseif (delta1 >0)&&(delta2 >0)% 历史最佳值<领域最佳值>当前值
        tubu(now_tabu_length,:) = var_value;% 将当前值写入禁忌表
        now_tabu_length = now_tabu_length + 1;
        %将领域最佳值更新为当前值
        func_value = best_near_func_value;
        var_value = best_near_var;
        %将领域最佳值更新为历史最佳值
        best_func_value = best_near_func_value;
        best_var = best_near_var;
        best_iteration = iteration;
    elseif (delta1 >0)&&(delta2 <=0)% 历史最佳值>领域最佳值>当前值
        not_in_tabu = true;% 不在禁忌表标志
        for i = 1:tabu_length
            for j = 1:var_num% 邻域最佳变量对比遍历禁忌表所有变量
                if abs(best_near_var(j) - tabu(i,j)) <= eps
                    in_tabu(j) = 1;
                else
                    in_tabu(j) = 0;
                end
            end
            if all(in_tabu)% 当前一组别变量中存在邻域最佳变量
               not_in_tabu = false;
               break;
            end
        end
        % 不在禁忌表标志
        if not_in_tabu
            % 将当前值写入禁忌表
            tubu(now_tabu_length,:) = var_value;
            now_tabu_length = now_tabu_length + 1;
            %将领域最佳值更新为当前值
            func_value = best_near_func_value;
            var_value = best_near_var;
        end
    end
    if now_tabu_length>tabu_length
        now_tabu_length = 1;% 防止禁忌表溢出
    end    
    all_best_func_value(iteration) = best_func_value;% 记录历史最佳值
end
%% 处理结果
% 打印结果
max_x = best_var(1);
max_y = best_var(2);
max_ans = target(max_x,max_y);
disp 最优迭代次数
disp(best_iteration)
disp 最大值位置
fprintf("x=%f\ny=%f\nz=%f\n",max_x,max_y,max_ans);
% 可视化结果
figure('position',[200 200 1280 480])
subplot(1,2,1);hold on;
[x,y]=meshgrid(-4:0.1:4);
z=target(x,y);
mesh(x,y,z);
title('Target Function')
xlabel('x');ylabel('y');zlabel('z');
view(-37.5,30)
pause(0.1)
scatter3(max_x,max_y,max_ans,'r*');
subplot(1,2,2);hold on;
x = 1:max_iteration;
y1= all_best_func_value;
try
    assert(length(x)==length(y1))
    plot(x,y1,'-','LineWidth',1.2)
    xlabel('x');ylabel('y');
    title('迭代数据')
catch
    disp("plot failed")
    disp(size(x))
    disp(size(y1))
end